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ABSTRACT 

It is widely believed that star clusters form with low star formation efficiencies. With the onset 
of stellar winds by massive stars or finally when the first super nova blows off, the residual 
gas is driven out of the embedded star cluster. Due to this fact a large amount, if not all, of 
the star s become unbo und and disperse in the gravitational potential of the galaxy. In this 
context. Ikroupal d2002l) suggested a new mechanism for the emergence of thickened Galactic 
discs. Massive star clusters add kinematically hot components to the galactic field populations, 
building up in this way, the Galactic thick disc as well. In this work we perform, for the first 
time, numerical simulations to investigate this scenario for the formation of the galactic discs 
of the Milky Way. We find that a significant kinematically hot population of stars may be 
injected into the disk of a galaxy such that a thick disk emerges. For the MW the star clusters 
that formed the thick disk must have had masses of about 10 6 M Q . 

Key words: Galaxy: disc — Galaxy: formation — Galaxy: kinematics and dynamics — 
galaxies: star clusters: general — methods: numerical 



1 INTRODUCTION 

The formation of the Milky Way (MW) galaxy is a mystery un- 
solved yet. Different models are trying to explain what were the 
initial conditions that lead to the actual structure of the MW. It is 
commonly accepted that the structure of the Milky Way, and other 
comparable disc galaxies, can be divided into three main compo- 
nents, the bulge, the galactic spheroid and the disc. The central 
bulge has a mass of « 10 10 M0 and a characteristic radius of about 
1 kpc. The galactic spheroid, which is also called the stellar halo, 
has a mass of w 3.7± 1.2 X 10 8 M teell et alj2008h and its mass 
is mostly confined within the the Solar Radius. The stellar halo con- 
tains globular clusters, dwarf satellites and their tidal streams, that 
add up to a mass of about 10 6 ~ 7 Mq. 

For the MW, the disc can be subdivided into at least two parts: 
the thin disc with a mass of about Mdisc = 5 x 10 10 Mq, that 
has exponential radial a nd vertical scale leng ths of approximately 
h R = 2.3 1 ± 0.6 kpc faammer et all 120071) and h z « 300 pc 
jjuric et all ,2008) respectively. The other part is the thick disk 
which has scale l engths of fethd.R = 4.1 ± 0.4 kpc and /ithd.z = 
0.75 ± 0.07 kpc dde Jong et ai]|20ld) . Near the Sun, the thick disc 
comprises about 6 per cent of the thin disc mass, so that the thick 
disc mass amounts to Mthd ~ 0.2 — 0.3 x Mdisc- The thick 
disc is made up mostly of low-metalicity ([Fe/H] ^ —0.4) stars 
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that have a velocity dispersion perpendicular to the disc plane of 
o"z,obs ~ 40 pc Myr~\ compared to the significantly smaller a z of 
the th in disc, which vari es from about 2 — 5 pc Myr~ 1 for the young 
stars dFuchs et"al I. I2OO1I) . These authors measured the velocity dis- 
persion in the solar neighbourhood as a function of age of the stars. 
The oldest stars in their sample (CNS4) have about 25 pcMyr -1 
for 10 Gyr old stars. But then this value might be 'contaminated' 
by thick disc stars. 

Several mechanisms have been proposed to explain the for- 
mation of th e thick disc in gala xies. One of these mechanisms was 
proposed bv lAbadi et al. who suggest that the formation of 

the thick disc is the direct accretion of stars from disrupted satel- 
lites. The process of accretion occurs appro ximately at coplanar 
orbi ts. Another explanation was suggested by iRoskar et al.l J2008h 
and lSchonrich & Binnevl j2009t) . who consider the process of ra- 
dial migration of the stars. In this mechanism, the stars which end 
up in the thick disc are trapped onto a resonant co-rotation with 
spiral arms and may migrate inwards and outwards along the spi- 
ral waves. This process conserves angular momentum and does not 
lead to significa nt heating of t he dis c . Another possible scenari o is 
proposed by e.g. Ouinn et all dl993l) . iKazantzidis et alj j2008l) and 
iVillalobos & HelmT i 20081 ) and consists of the thickening of a pre- 
existing thin disc through minor mergers. The thick disc is formed 
by the dynamical heating that is induced by satellites merging with 
a prim ordi al, rotationally s uppor ted thin disc. Finally, iBrook et all 
and lBornaud et all d2007[) suggest that the formation of the 
thick disc is triggered in situ. The process of star formation occurs 
during/after gas rich mergers. 

Each model explains different aspects and has its own impli- 
cations for the disc's kinematical and chemical properties. In the 
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scenario of lAbadi et al.l d2003h . it is possible to distinguish two dif- 
ferent dynamics for the orbits of the stars, and they allow one to 
identify two different components of the disc. The thin disc is a 
kinematically cold component with stars on circular orbits, and the 
thick disc is a kinematically hot component with stars with orbital 
parameters transitional between the thin disc and the spheroid. The 
nature of the population of stars in the thick disc comes from the 
tidal debris of satellites, whose orbital plane is roughly coincident 
with the disc. Their orbits were circularized by dynamical friction 
prior to their complete disruption. 

However, there is evidence for enhanced a-elements abun- 
dances in thick-disc stars. This indicates a short star formation 
time-scale in which enrichment is dominated by Type II sup ernovae 
(SNe II) jAlves-Brito eTail |20ld : iKobavashi et all 1201 ll) . rather 
than the slower time-scale expected for dwarf galaxies. In this case, 
one needs to consider the existence of an active epoch of gas-rich 
mergers in the p ast history of galax y, with the bulk of the thick disc 
forming in situ ( iBrook et al.ll2005l) , and stron g stellar scattering of 
clumps formed by gravitational instabilities dBornaud et alj|2007l) 
(rather than being accreted from satellites) to explain these obser- 
vations. 

During the last years understanding of star formation and star 
cluster formation has progressed. Now we know tha t almost all 
stars form in a clustered mode (e.g. lLada & Ladal2003l and follow- 
up publications). The star clusters form with low star formation ef- 
ficiencies and, therefore, they lose a large part of their stars that ex- 
pand outwards when resi dual gas is exp elled by the action of mas- 
sive stars. In this context, iKroupal J200% suggested another mecha- 
nism for the formation of the thick disc. Considering the star cluster 
formation process described above, massive star clusters may add 
kinematically hot components to galactic field populations, build- 
ing up, in this way, bot h galactic discs . 

According to the IKroupal d2002h model, the velocity disper- 
sion of the stars in expansion (the new field population) is related 
with the mass of the initial star cluster, wi th its radius an d also with 
the efficiency of star formation (Eq. 1 in lKroupall2002l) . It is pro- 
portional to the square root of the mass of the star cluster, such that 
heavier star clusters will generate a field population that presents 
higher dispersion velocities. Thus, a simple estimate suggests that 
clusters of mass near 10 M@ may contribute new field popula- 
tions that have a velocity dispersion of about 40 pc Myr -1 which is 
similar to the velocity dispersion of the thick disk. 

In this work we perform numerical simulations to investigate 
this scenario for the formation of both galactic discs of the Milky 
Way. To accomplish this investigation we use star clusters of dif- 
ferent masses and we consider their disruption under distinct low 
star formation efficiencies (we refer to such clusters as popping star 
clusters). We place the star cluster on an orbit around the Galactic 
Centre at the Solar Radius of 8.5 kpc and observe how the stars of 
these star clusters distribute themselves in the potential of the MW. 

In the next Section we describe the detailed setup of our in- 
vestigation and describe our results in Section [3] We discuss our 
findings in the last Section {4). 



2 SIMULATIONS 
2.1 Code Description 

We simulate the popping star cluster using the particle mesh 
code SUPERBOX dFellhauer et al.ll2000T) . which has moving high- 
resolution sub-grids, which stay focused on the star cluster. These 



sub-grids provide high spatial resolution at the place of interest. A 
particle-mesh code neglects by default close encounters between 
the particles (which are rather representations of the phase-space 
than actual single stars) and is therefore called collision-less. With 
a collision-less code it is possible to simulate galaxies without hav- 
ing to use the actual number of stars (10 10 ). A code as SUPERB OX 
is in principle unsuitable to simulate star clusters, as with them two- 
body relaxation effects are important and govern their evolution. 
In our study we dissolve the star clusters immediately, dispersing 
their stars within a short time (a few crossing times) into the over- 
all potential of the MW, therewith being in a regime where near 
encounters between stars do not play any role. Now the ability of 
SUPERBOX to model particles with arbitrary masses gives us the 
possibility to sample the phase space (e.g. by using more particles 
than actual stars in the star cluster) more precisely. 

The code is fast and resource-efficient, i.e. it requires a small 
amount of computer memory. Therefore, this code enables its user 
to simulate objects with millions of particles, with high resolution, 
and it can run on normal desktop computers. 

SUPERBOX has two levels of high resolution sub-grids. The 
highest resolution grid has a resolution (i.e. cell-length) of 0.2 pc 
and covers the initially dense SC completely. The medium resolu- 
tion grid has a cell-length of 40 pc and covers the complete z-height 
of interest. Finally the outermost grid covers the complete orbit of 
the SCs around the centre of the MW with a resolution of 0.4 kpc. 

Regarding the fixed time-step of SUPERBOX we choose ini- 
tially a value to ensure that the SCs are kept stable without further 
influence. Therefore the time-steps were different for the different 
masses ranging from 0.03 Myr for the 10 M0-clusters, via 0.01 
and 0.006 Myr to 100 yr (0.0001 Myr) for the most massive clus- 
ters (10 7 Mo). As these tiny time-steps require very long simula- 
tion times until we finally would reach 10 Gyr, we check the simu- 
lations and stop them as soon as the SCs are completely dissolved. 
When all stars are unbound and just traveling under the influence of 
the global potential, we restart the simulations with a larger time- 
step of 0.5 Myr. 



2.2 Set-up 

We consider the following set-up for our simulation. The potential 
of the MW is modeled as an analytical background potential con- 
sisting of a Hernquist sphere to generate the bulge, 
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using M b = 3.4 x 10 10 M Q and a = 0.7 kpc, the Miamoto-Nagai 
model to mimic the disc, 
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with M d = 10 M s , b = 6.5 kpc and c = 0.26 kpc. Finally, we 
use a logarithmic potential to account for the rotation curve of the 
MW disk, 



<£>haio(r) = y m(r 2 



(3) 



with vo = 186 km s~ x and d — 12 kpc. The superposition of these 
components gives a good analytical representation of the Milky 
Way potential today. The potential of the MW is kept constant 
throughout the simulation. This is highly idealised as we expect 
the Galaxy to grow and evolve over this period of time. We discuss 
this issue further in Sect.|4] 
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Table 1. Initial conditions of the Plummer spheres which represent the star 
clusters for each simulations. The first column shows the initial Plummer 
radius of the star clusters, the second the initial mass of the sphere and the 
third the initial characteristic velocity dispersion. Finally the fourth and fifth 
columns represent the SFE and the final mass after gas expulsion. 
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The single cluster is represented by a |Plummer1(ll91 ll) sphere. 
Even though, this might not be the most ideal representation of 
a young embedded star cluster, the Plummer distribution has the 
advantage that all quantities are easily accessible with analytical 
formulas. And it is i n fact a better representation of a young, newly 
formed star cluster dKroupa|[2008l) . which has not yet been tidally 
truncated, than e.g. a King profile. The Plummer profile has two 
free parameters to choose, namely the Plummer radius, which we 
keep fixed at R p \ — 1 pc, and the total mass of the Plummer sphere 
- a parameter which is varied in this study. All other quantities (e.g. 
crossing-time, velocity dispersions, distribution of velocities) are 
then simple functions of the two input parameters. The choice of the 
small Plummer r adius is based on observations of embedded star 
forming regions jKroupall2008h . Furthermore we choose a cutoff 
radius of 5 pc. The choice of the cut-off radius is justified by the fact 
that beyond 5R P \ there is only less than 2 per cent of the total mass 
of the Plummer model missing. Each cluster has an initial mass of 
M cc i = 2.2 • 10 4 , 2.2 • 10 s , 2.2 ■ 10 6 and 2.2 • 10 7 M respectively 
leading to an initial crossing time of 0.6 Myr, 0.2 Myr, 0.06 Myr 
and 0.001 Myr and is represented with 2, 500, 000 particles. 

This configuration, placed on a circular orbit at the Solar Ra- 
dius of 8 kpc, suffers from the process of gas expulsion. This pro- 
cess is mimicked by artificially reducing the mass of each particle 
according to the SFE used. We investigate a range of possible SFEs 
namely 0.4, 0.2 and 0.01. 

Finally, we also use two distinct ways, how this mass-loss pro- 
gresses. We perform simulations with an instant mass-loss and sim- 
ulations where the mass in gas is lost during the time-interval of one 
crossing time of the star cluster. 

As the simulation time we choose a time-span of 10 Gyr to 
ensure that the particles have enough time to disperse into a stable 
configuration in the Galactic potential. 



3 RESULTS 

3.1 2 -distribution of the stars 

In Fig.[T] we show the number of stars with respect to their z-height 
after 10 Gyr of simulation time. Therefore, we count all stars in a 



Table 2. Results of the simulations. The first three columns represent the 
initial mass of the star clusters, the SFE used and finally if the gas was 
expelled instantaneously or over a crossing-time 1 (GET = gas expulsion 
time). The fourth and fifth column represent the scale height of the thin and 
thick disc, respectively, after 10 Gyr of simulation. The last line represents 
our combined model (see main text). 
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vertical cylinder of 40 pc radius around the Sun. To increase the 
statistical significance and also because we do not know the ac- 
tual position of the star cluster with respect to the Sun (except that 
they have by design the same Galacto-centric distance), we actu- 
ally count the stars in an 80 pc wide (7.96 «; R = y/x 2 + y 2 ^ 
8.04 kpc) ring around the galaxy and plot the z -height of the stars 
in Fig. Q] We choose to simulate until 10 Gyr to ensure that the 
objects had enough time to evolve slowly into an equilibrium. The 
SCs at this stage are completely dissolved in almost all of the cases. 

The different rows of Fig. [T]correspond to the three different 
SFEs we adopt. The top row shows the results of our SCs having 
40 per cent of SFE, the middle row the results of the 20 per cent 
simulations and finally the bottom row the results obtained with the 
extremely low efficiency of only 1 per cent. The left column shows 
simulations where the gas-removal time was one crossing-time of 
the star cluster, while the right column represents the final distribu- 
tion of stars when the gas was removed instantaneously. The curves, 
in each panel, from left to right correspond to the star distributions 
of star clusters which had an initial mass of M ec i = 2.2 x 10 4 M 
(black on-line), 2.2 x 10 5 M Q (red on-line), 2.2 x 10 6 Mq (green 
on-line) and finally 2.2 x 10 7 M Q (blue on-line). 

First we check our results for a dependency on the under-lying 
analytic disc potential (Miamoto-Nagai model). For this reason we 
compare different fitting functions to our distributions. One would 
expect that if the under-lying potential had any influence it would 
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Figure 1. The ^-distribution of particles of the dissolved star clusters after 10 Gyr of simulation time. From top to bottom we show the results of the different 
SFEs (0.4, 0.2 and 0.01). Left panels are for the gas-expulsion during a crossing-time and right panels show the results for instantaneous gas-expulsion. Inside 
each panel the curves from left to right represent the results of the different initial masses of the star clusters namely M ec \ = 2.2 X 10 4 Mq (circles black 
on-line), 2.2 X 10 s Mq (squares red on-line), 2.2 X 10 6 Mq (triangles green) and finally 2.2 X 10 7 Mq (crosses blue). 



be visible in form of a similar distribution of the dispersed stars, 
i.e. in our case a Plummer like distribution, as the Miamoto-Nagai 
potential is a flattened version of the Plummer potential. But our 
results showed that the best fitting functions were either a single 
or a double exponential. Therefore we are quite confident, that our 



results are not influenced by the particular choice of the analytical 
potential. 

The results show no difference with respect to the method of 
gas- (mass-)loss we applied. The final configurations are similar 
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Figure 2. Distribution of stars in the .z -direction if we add our different mass 
clusters according to a standard ICMF. Fitting line is a double exponential. 



irrespective if we remove the gas instantaneously or more smoothly 
during a crossing-time of the star cluster. 

But we see immediately a strong dependency on the initial 
mass of the star cluster. Low mass clusters (2.2 x 10 4 M Q and 
2.2 x 10 5 M©) disperse their stars only close to the Galactic plane 
and therefore can only contribute to the build-up of the thin disc. 
This changes as soon as we look at the high mass clusters (2.2 x 
10 6 M and 2.2 x 10 7 M Q ). These clusters are able to distribute 
their stars out to 1-4 kpc in z-height and therefore could be the 
building blocks of the thick disc. 

Finally the SFE has only a second order effect on the final 
distribution. The lower the SFE the wider is the spread of the fi- 
nal particle distribution. The reason is quite simple as a lower SFE 
means that the initially binding potential of the embedded cluster 
is deeper such that the stars, which amount to a mass M say, move 
more rapidly within it. The same stellar mass, M, but confined 
by a less-massive gas potential (being equivalent to a higher star- 
formation efficiency), has a smaller velocity dispersion before the 
gas is removed. One also has to note that even though a SFE of 0.4 
should in principal lead to a surviving remnant star cluster, for all 
masses except our highest mass models, these remnants were not 
massive enough to survive the whole 10 Gyr of simulation time. 
With a SFE of 0.4 the cluster still 'pops' but instead of dissolv- 
ing comple tely a small fraction of stars re-virialises into a bound 
star cluster. iKroupa et alj d200ll) performed high-precision Nbody 
modelling of this process explaining that a pre-Orion Nebula Clus- 
ter evolves into a Pleiades type cluster whereby 2/3rds of the pop- 
ulation "pop away". The majority of stars thus gets distributed as 
unbound stars inside the MW potential. The fraction of stars which 
remains bound is a function of the SFE. For the low-mass cases 
this bound remnant does not survive the tidal forces of the MW for 
very long. Only the more massive remnants stemming from mas- 
sive initial clusters are able to survive the full 10 Gyr. The future 
fate of these remnants is of no concerns in our study, as our code is 
not able to describe the internal dynamics of star clusters over long 
time intervals correctly. 

As stated above we then fitted exponential functions to the 
obtained star distributions in ^-direction. For all models we tried to 
fit a single exponential of the form: 



Table 3. ly-velocity dispersions of our dissolving star clusters. The first 
three columns are the same as in Tab.|2] We then give the fitted values for the 
dispersion using the method described in the main text for the thin and thick 
disc component. The last column is the FWHM of the whole distribution. 
Low mass star clusters, which only contribute to the thin disc have only one 
value fitted. Last line is the result of our combined model using an ICMF. 
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p(z) = p ,thin • exp , (4) 

where ft z ,thin represents the z-scale height of the thin disc and com- 
pared the results with the ones obtained if we used a double expo- 
nential function of the form: 

p(z) = po.thin ■ exp ( — — — I + po.thick • exp ( — — — I , (5) 

\ftz,thin / \1z, thick / 

where /i Zl thick is the z-scale height of the thick disc. The fitted val- 
ues of /i z ,thin and /i z , thick are shown in Tab. [2] For the low mass 
clusters the best fit was always the single exponential fit, while for 
the high mass clusters the distribution was better fitted with a dou- 
ble exponential. In the table we only give the values for the best 
fit of the two fitting functions adopted. Again it is clear that the 
method of removing the gas has no influence on the final distri- 
bution. But we see a clear trend towards larger scale-lengths by 
increasing the cluster mass and by lowering the SFE. For an exam- 
ple, in the case M ec i = 2.2 x 10 6 M© and a SFEs of 40% or 20%, 
the z scale-height of the thin disc is about 14 pc. Using a SFE of 
only 1 per cent gives a z-scale height of 16-18.2 pc. The same type 
of behavior can be observed for the thick discs. 

We should make clear that the actual distribution of stars is 
neither a perfect single nor a perfect double exponential distribu- 
tion. With the low-mass clusters we only see a thin distribution 
like a 'peak' and only few stars in a sort of 'envelope' but with 



6 Assmann et al. 



low z -heights. Therefore a single exponential fits the data well. For 
high-mass clusters we see this 'peak' as well together with a very 
extended structure (the 'envelope') reaching out to large z-heights. 
These clusters are able to spread out their stars to large radii. For 
those distributions a double exponential fits the data better. 

Even though none of the values from Tab. [2] come anywhere 
close to the values of the thin and thick disc of the MW, it is as- 
tonishing that our high-mass clusters distribute their particles auto- 
matically into a distribution which is best fitted by a double expo- 
nential, i.e. show the same shape as the actual distribution of MW 
stars. The low numbers can easily be understood by the fact that 
the MW is not made out of one single dispersed star cluster at one 
certain time and on one certain orbit. We rather have an overlay of 
star clusters forming with different masses at different times and on 
all radial distances. Furthermore, we do not take any further effect 
into account which could enhance the vertical spread of the star 
particles, e.g. scattering with giant molecular clouds or spiral arms, 
adiabatic heating due to the growth in mass of the MW disk. 

3.1.1 Combined model 

To at least asses the influence a generation of star clusters would 
have onto the results we build our combined model by assum- 
ing an initial cluster mass function (ICMF). The ICMF is a fun- 
damental property of the process of star formation in galaxies 
dKroupa & Boily||2002h . It gives the number of star clusters in a 
certain mass interval dAl cc \ which form during one event of star 
formation, 



N(M e , 



oc M oc fdM ec i, 



(6) 



where /3 is the spectral index of cloud mass spectrum. In this work 
we use P = 2, because this value was measur ed for young star 
clusters in the Antennae Jwtiitmore et al.lll999h and generally by 
IWeidner & KrouptJ Jiooj) . Because we use a particle-mesh code, 
where particles represent phase-space elements and not single stars, 
all of our model clusters have the same number of particles, irre- 
spective of their mass. We therefore are able mimic an ICMF with 
a power law index of p — 2 by co-adding our results using a power 
of P — 1 = 1 only, to ensure that every particle has the correct 
weight in the combined model. 

The density distribution of the disrupted star clusters at 
10 Gyr using the ICMF is shown in Fig [2] In this figure we con- 
sider star clusters with a SFE of 20% and where the gas expulsion 
happens during a crossing time. One clearly sees that the shape 
of the density distribution has a similar resemblance as the thin and 
thick disc of the MW. We fit a double exponential function given by 
Eq.|5]and obtain values for the z scale-heights of /i z ,thin = 50 pc 
for the thin disc and ft z , thick = 500 for the thick disc (see also last 
line of Tab.[2]l. These values are now much closer to the real values 
of the MW and therefore we believe that by adding also additional 
scattering and multiple star formation events, we may be able to 
account for the real discs of the MW. 



3.2 VK-velocity distributions 

In addition to studying the shape of the z-distribution of the stars, 
we can also use our simulations to look at the velocity space. Ta- 
ble[3]shows the velocity dispersion obtained for all our simulations 
by fitting straight lines to the probabi lity plot of W -velocitie s (ve- 
locity in z-direction), as described in iBochanski et~aH d2007l) . The 
probability plot is a graphical technique where we plot the cumu- 
lative probability distribution in units of the standard deviation of 
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Figure 3. In the top panel we show the W-velocity distribution of one of 
our high mass models (see main text). One clearly sees the non-Gaussianity 
of the distribution. The lower panel shows the probability plot for this distri- 
bution with two straight lines fitted to the inner and outer parts respectively, 
giving the velocity dispersion for the two parts. 



the distribution. When the distribution is a Gaussian, it shows as a 
straight line in the probability plot with a slope corresponding to 
the standard deviation. The t/-axis-intercept is equal to the median 
of the distribution. When the distribution is non-Gaussian, the plot- 
ting line deviates from a straight line. We are therefore able to fit 
two separate lines to the inner and outer parts of the distribution, 
resembling the velocity dispersion for the thin and thick disc com- 
ponent of our models. Again, we note that our low-mass models 
only contribute to the thin disc component. 

As example we show the velocity distribution and its probabil- 
ity plot for a high mass model in Fig. [3] It shows a star cluster with 
a SFE of 20 % and a mass of 2.2 x 10 6 M . Clearly, the velocity 
distribution of the stars is not a Gaussian at all. Also two Gaussians 
will not represent this kind of distribution fully. But it shows clear 
similarities with the ve locity distribution of stars found in the MW 
jBochanski et alj|2007l) . So even though we have an overlay of an 
undefined number of Gaussians we perform the same analysis as 
observers and fit straight lines to the inner and outer parts of the 
probability plot. The probability plot is well fitted with two lines: 
low dispersion velocities (<| la ), from the population of stars 
confined to the thin disc, and high dispersion velocities (>| la ), 
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Table 4. Ratio between initial velocity dispersion of the star clusters and 
the final velocity disspersion of the dispersed stars. Shown are the models 
with a SFE of 0.01 and instanteneous gas-expulsion only. 



initial mass 


final <7 


initial a 


ratio 


[M ] 


[kms- 1 ] 


[kms" 1 ] 




2.2 X 10 4 


2.0 


5.3 


0.38 


2.2 X 10 5 


2.8 


16.7 


0.17 


2.2 X 10 6 


27.3 


52.8 


0.52 


2.2 X 10 7 


69.8 


166.8 


0.42 



3.2.1 Combined model 

As with the z— distribution we now investigate how our results 
change if we adopt an ICMF. We combine our models the same 
way as done for the ^-distribution and show the combined velocity 
distribution in the top panel of Fig. [4] (and also in the last line of 
Tab. [3j. The derived probability plot to this distribution is shown 
in the lower panel of Fig. [4] Fitting straight lines to the inner and 
outer parts give velocity dispersions of 1.3 and 80 kms -1 for the 
thin and thick disc component respectively. 

As we can see, the velocity dispersion is much larger in the 
thick disc than in the thin disc. Comparing o ur values of the so- 
lar neighborhood with published observ ations (iBensbv et alj2003t 
IVallenari et"Zll200rj : fVeltz et al]|2008h . which report <r w ,thick = 
38 kms" 1 and <r w ,thin = 16 kms -1 for the thick disc and thin 
disc respectively, shows that we underestimate the thin disc values 
quite substantially and overestimate the value for the thick disc. An 
explanation could be that we do not take any other scattering mech- 
anisms into account, which would elevate our thin disc dispersion. 
For the thick disc it could probably be, that there never were any 
SCs with 10' Mq formed in the history of the MW. These high- 
mass clusters are mainly responsible for the elevated velocity dis- 
persion in the thick disc component of our models. 



Figure 4. Same as Fig.[3]but now for our combined model using an ICMF 
as described in the main text. 



where the population of the stars corresponds to the thick disc. In 
our example, the slope for the velocity dispersion for the stars of the 
thin disc is 4.01 km s _1 and for the thick disc stars 20.5 km s~ . 

Additionally we measure the full width half maximum of the 
velocity distribution and give these values as the last entry in Tab.[3] 

If we compare the velocity dispersions of the initial star clus- 
ters with the dispersions we measure of the dispersed stars we see 
the following effects: 

• The dispersion of the dispersed stars is always significantly 
lower than the initial dispersion inside the star cluster. 

• The ratio of final to initial velocity dispersion decreases with 
increasing mass (i.e. initial dispersion) of the star cluster. 

• This fraction is higher for the stars in the thick disc compo- 
nent. 

• We get higher dispersions for a given initial cluster mass, if 
the SFE is lower. The reason is that the stars at lower SFEs have a 
smaller potential barrier to cross. 

• If the initial mass of the star cluster is about 10 6 Mq the ve- 
locity distributions are better fitted by two Gaussians than just one; 
similar to the z-distributions. 



3.3 Of Christmas trees and azimu thai fish 

One interesting result can be observed in our simulation when we 
focus on the z-distribution of the stars at early stages (< 1 Gyr) 
of the SC evolution. For an example, in Fig [5] we show the z- 
distribution of the stars, after 270 Myr, for the simulation with an 
initial mass of 2.2 X 10 7 M and a SFE of 0.2. 

The shape of the z-distribution looks like a Christmas tree 
(this effect was seen in our simulations shortly before Christmas 
2009). The origin of this shape can be explained by considering 
the dynamical processes during the SC evolution. The star cluster 
is orbiting on a circular orbit around the Galactic Centre. After gas 
expulsion, some of the stars acquire higher expanding velocities 
and form the thick disc. We show in Fig|6]in the top left panel the 
W-velocity distribution as a function of the distance between 8 and 
9 kpc from the Galactic Centre after 270 Myr (i.e. about one orbital 
period). The velo cities show di s crete l ayers. A similar behaviour is 
also described in Kii pper et al.l ( 120 101) where the tidal tails of SCs 
show density enhancements with regular spacings stemming from 
the turn-around points of the escaped stars on epicycles. We sus- 
pect that our Christmas tree has a similar explanation and is due to 
the turn-around poin ts of stars on their z-oscillations (z-epicycle). 
I Kiipper et alj d2010h report that their star clusters have to be suffi- 
ciently long lived to show this behaviour in their tidal tails, while 
we see this effect very early in the evolution of our models. The 
reason for th i s diff erence is quite simple. In the simulations of 
iKiipper et all ( l2010l) the star clusters are in virial equilibrium and 
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Figure 5. Example for a Christmas tree as found in all of our simulations. 
It is a peculiar transient feature seen only after one revolution around the 
galaxy RS 270 Myr (also the simulation time shown in the figure). It disap- 
pears after a few revolutions more. 



loose stars through tidal and two-body effects quite slowly. There- 
fore i t needs a long time fo r the tidal tails to form and evolve. Thus, 
in the lKtipper et alj j201Ch case the epicyclic motions of stars evap- 
orating from star clusters are nearly coherent or in-phase since the 
stars leave the cluster with very similar and small velocities. The 
popping clusters, on the other hand, expel stars with a larger range 
of velocities within a short time-span. 

The middle and lower left panel of Fig. [6] show that we only 
see this phenomenon in the W^-velocities, while U (radial veloc- 
ity) and V (azimuthal velocity) show a simple bulk motion (except 
for a few stars which lag behind). This special appearance of the 
Christmas tree gets later erased and transformed into a smooth dis- 
tribution as described in Sect. 13. II Also the velocities in W do not 
show the quantized behaviour any longer as we can see in the top 
right panel of Fig. [6] The transient tree feature more or less disap- 
pears when the stars which form the leading arm of the tidal tails 
wrap all around the whole orbit and start to overlap with the trailing 
stars filling the gaps in velocity space. A deeper theoretical expla- 
nation for this astonishing feature is not part of this paper and will 
be dealt with in a future publication. 

In the right panels of Fig. [6] we show the particle velocities 
after 10 Gyr of simulations. While we clearly see that the quan- 
tized nature in the V^-component has washed out (see above) a 
new feature has appeared. In all our simulations we see over- and 
under-densities in the azimuthal V-velocities. This quantization is 
not as sharp as formerly seen in W but still clearly visible. We 
call this feature azimuthal fish, because of its form. We believe that 
this quantization is the direct counterpa rt of the density enh ance- 
ments in tidal tails which are reported bv lKupperetailfeoiOl) . The 
only difference between their and our models is that they inves- 
tigate slowly evolving globular clusters (GCs), which loose their 
mass slowly over long time-intervals, while our models have dis- 
solved rapidly and the stars spread out much faster and wider than 
seen in lKtipper et alj {201Q). We will investigate this problem fur- 
ther and report about it in a follow up paper. 



We compute models of dissolving embedded star clusters on circu- 
lar orbits around the MW. The aim of our study is to investigate if 
SCs of different mass, SFE and gas-expulsion time can support not 
only the distribution of stars and their velocities in the thin disc of 
the MW but also in the thick disc. 

We also co-added the results of our simulations to show the 
effect of a whole population of SCs, following a standard ICMF, to 
the distribution of stars in the discs of the MW. 

We can summarize the results of our models in the following 
points: 

• The stars of dissolved SCs form density distributions in the 
z-direction which can be best fitted by exponential profiles. 

• While low-mass SCs (masses of the order of 10 4 and 10 5 Mq) 
contribute only to a component similar to the thin disc of the MW, 
high mass clusters(masses > 10 Mq) show distributions which 
can be described best by two exponential profiles and therefore 
their stars contribute into two components similar to the thin and 
thick disc of the MW. 

• The velocity dispersion of the stars distributed into the thin 
component show a much lower velocity dispersion than the actual 
dispersion of the thin disc of the MW. This may be due to the fact 
that we do not take any other mechanisms to enhance the velocity 
dispersion into account (e.g. spiral arms, giant molecular clouds, 
other SCs, ...) 

• The thick component of our model shows a velocity dispersion 
which is higher than observed in the MW thick disc. We therefore 
conclude that inside the MW disc SCs with masses comparable to 
10 7 Mq might never have formed. The thick component of our 
10 6 Mq model has enough velocity dispersion to explain t he thick 
disc. T his is in nice agreement with the analytical results of lKroupi] 
j2002l) . 

We need to discuss these findings a bit further. On the first 
glance it seems to be rather odd that stars of one dissolved star clus- 
ter should spread into two distinct distributions with two distinct 
velocity dispersions. The density distribution of all our dissolved 
clusters exhibits a very distinct profile showing a 'peak' around 
small z-distances and a kind of 'envelope' of large z-distances. 
The transition between these two parts is rather continuous and not 
sharp. We see that for low-mass clusters the 'peak' is more pro- 
nounced than the 'envelope' and therefore they can be fitted best 
by a single exponential profile. Also the stars in the 'envelope' do 
not extend to z -heights comparable to the thick disc. For high-mass 
clusters we see a very substantial 'envelope', reaching out to high 
z-distances, and therefore are better described by two exponentials. 
The same is true for the distribution of M^-velocities. The distribu- 
tion is similar in shape to the distribution of velocities found in 
stars in the solar neighbourhood. It is neither a single nor a double 
Gaussian but rather a continuous overlay of many Gaussians. Nev- 
ertheless, the usual observational procedure to analyse these distri- 
butions is to fit two Gaussians to the probability plot and assigning 
two velocity dispersions to it - one for the thin disc component and 
one for the thick disc component. So we followed the same method. 
Again for low-mass clusters there are no stars with high velocities 
and therefore we only assign one velocity dispersions to their dis- 
tributions. 

Another issue, which needs to be discussed, is the constant 
analytical potential used in this study. It is clear that in the distant 
past the disc(s) of the MW were much smaller in mass and maybe in 
size. A 'popping' star cluster would spread its stars to much larger 
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Figure 6. Evolution of a model SC in velocity space. Left panels show the velocity distributions after 270 Myr (rj one revolution round the Galaxy). Right 
panels show the distributions at the end of the simulation (10 Gyr). From top to bottom we show W against R, W against U and finally W against V. Top left 
panel shows clearly the quantized nature of the W-velocities which lead to the so-called Christmas-tree effect. At the lower right panel we see the 'azimuthal 
fish'. 



z-heights much more easily. So it would rather help our scenario 
to explain the thick disc. Furthermore the adiabatic heating of stars 
due to the growth in mass of the MW will enhance the velocity 
dispersion of the stars and bring them onto orbits with higher z- 
distances. 



Besides, it is well known that the thick disc and thin disc 
components have a different chemical history. The ratios of a- 
(0,Mg,Si,Ca, and Ti) elements, for example, are higher for stars 
in t he thick disc than f o r stars of the thin disc at given metallic- 
ity teensbv et aT]|2007l : iKobavashi et all 1201 lh . This indicates a 
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short star formation time-scale in which enrichment is do minated 
by Type II supernovae (SNe II) l lAlves-Brito et al.ll2010l) . Super- 
novae events guide the gas expulsion process, that together with 
all the other scattering mechanisms, might also help to move stars 
from the 'peak' structure into the 'envelope' structure. This would 
explain the problem that we do not see stars with thick disc chem- 
ical properties amongst thin disc stars and, therefore, why we see 
different populations of stars with different chemical properties in 
the thin and thick disc of our MW now. Another possibility would 
be if the MW saw an early epoch of star formation, where the for- 
mation of low-mass clusters was suppressed (in comparison with 
the evidence for this from top-h eavy galaxy-wide in itial mass func- 
tions of stars in young galaxies TWeidner et al.l201oT) . 

We therefore are able to show that one does not need any 
merger scenarios to explain the structure of the MW disc. The sim- 
ple assumption that all stars form in SCs is enough to explain the 
thin as well as the thi ck disc o f our M W, thus confirming the ana- 
lytical calculat ions bv l Kroupal l2002h rather nicely. These models 
are applied by iKroupal 12003) as a possible explanation not only 
of the thick disk, but also of the observed but hitherto not under- 
stood secular heating of the thin disk. Here the notion raised is 
that the observed thickening of the MW thin disk with age may be 
due to a falling star-formation rate (SFR) since the past 10 Gyr, 
whereb y a small SFR implies th e formation of SCs with small 
masses dWeidner & Kroupall20o4) . 

While conducting our experiments we detected two new ef- 
fects within the distributions of the stars. After about one revolu- 
tion around the Galaxy the stars have spread in the z-direction in 
to a form that has striking similarities to a Christmas-tree. As the 
discovery was made shortly before Christmas we called the effect 
officially Christmas-tree distribution. At the same time the distri- 
bution of W-velocities seems to have only certain discrete values. 
This feature is transient and disappears after several revolutions 
about the MW. We suspect that this is the case when the leading 
and trailing arms wrap completely around the Galaxy and overlap 
at the former SC position. Then the velocities of the leading arm 
fall into the gaps of the velocity distribution of the trailing stars. 
The reason to see a peculiar distribution like the Christmas-tree in 
the first place is still not fully investigated and will be dealt with in 
a follow-up pap er. We suspect that w e see a similar behaviour in z 
as is reported by |Kupperetal] ( l2010l) along the tidal tails. 

The second effect showed up towards the end of our simula- 
tions. This effect consists of peculiar over- and under-densities in 
the V velocities which show a pattern that reminded us of a fish. 
We therefore called this effect the azimuthal fish. Even though the 
'tidal tails' of the dissolved star clusters are well mixed towards the 
end of the simulations so that we can not see any density variations 
in the positions, we believe that our fish pattern is directly related to 
the slowly evolving tidal tails of GCs and their density variations. 
As the quantized velocities associated with the Christmas-tree we 
see now the left-over velocity signatures of a similar effect in az- 
imuthal direction., while the positional counterpart is not visible 
due to the fast evolution of our models. Also this effect will be in- 
vestigated further and reported about in a following publication. 

As a final comment we state again that our models are able 
to reproduce the structural features of the MW. Our models add 
therefore another possible explanation of how the discs of the MW 
and similar galaxies may have acquired their properties. Which 
theory will be right - future observation might tell or even show 
that we might deal with a superposition of all theories at the same 
time. 
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